import numpy as np
import math
from pyXSteam.XSteam import XSteam
from tqdm import tqdm
import pandas as pd
import time
from scipy.io import savemat
import matplotlib.pyplot as plt
from matplotlib import rcParams
from iapws import IAPWS97
from pylab import mpl

plt.rcParams["font.family"] = ["SimSun", "Times New Roman"]  # 中文使用宋体，英文使用Times New Roman
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

"积分定义"
class Integrator:
    def __init__(self):
        pass
    @staticmethod
    def euler(dy, y, step_size, lower_bound=float('-inf'), upper_bound=float('inf')):
        y += dy * step_size
        y = np.clip(y, lower_bound, upper_bound)
        return y

def turbine_coff(D,P):
    k=P@np.linalg.pinv(D)
    return k

"汽轮机"
class steam_turbine:
    def __init__(self,ita):
        self.a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.integrator = Integrator()
        self.ita = ita

    "调节级"
    def regulating_level(self, Ts, Pst, Dst,k):
        """ 初始化"""
        # hst2=vst2=ht2=hst1=vst1=ht1=t=ht=hst=0
        """常数"""
        ptvwo = 14  # VWO工况下调节级后压力
        pstb2vwo = 16.67  # VWO工况下流经部分开启调节阀后的蒸汽压力
        Dstbvwo = 290.76  # 设计工况下流经部分开启调节阀后的蒸汽流量
        # ptd = 12.3  # 设计工况下调节级后压力
        # ptd = 12.3 #冬季
        # Dstd = 513.6  # 设计工况下主蒸汽流量
        fai = 1  # 喷嘴和动叶速度系数
        Kn = 1.3  # 蒸汽绝热指数
        Dstvwo = 581.5138889  # VWO工况下主蒸汽流量
        hst = self.a.h_pt(Pst*10, Ts)
        """调节级"""
        if Dst <= 290.7569444:
            Dt = Dst  # 调节级后蒸汽流量
            pt = Dst * k  # 变工况下调节级后压力
            pst2 = Dst * ((pstb2vwo - ptvwo) / Dstbvwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.a.h_pt(Pst * 10, Ts)
            vst2 = self.a.v_ph(pst2 * 10, hst2)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht = self.a.h_pt(pt * 10, tt2)
            t = self.a.t_ph(pt * 10, ht)
            # 调节级蒸汽温度tt
        elif 290.7569444 < Dst <= 436.1354167:
            """#部分开阀"""
            pt = Dst * k  # 变工况下调节级后压力
            pst2 = (Dst - 0.5 * Dstvwo) * ((pstb2vwo - ptvwo) / Dstbvwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.a.h_pt(Pst * 10, Ts)
            vst2 = self.a.v_ph(pst2 * 10, hst2)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht2 = self.a.h_pt(pt * 10, tt2)  # 通过部分开启调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值ht
            '''#全开阀'''
            hst1 = self.a.h_pt(Pst * 10, Ts)
            pst1 = 0.95 * Pst  # 蒸汽流经全开调节阀门后压力pst1
            vst1 = self.a.v_ph(pst1 * 10, hst1)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数
            vt1 = (pst1 / pt) ** (1 / nm) * vst1  # 通过全开调节阀蒸汽流经喷嘴动叶后比体积
            tt1 = (1000 * pt * vt1 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过全开调节阀蒸汽流经喷嘴动叶后温度
            ht1 = self.a.h_pt(pt * 10, tt1)  # 通过全开调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值
            """%混合"""
            Dt = Dst  # 调节级后蒸汽流量
            ht = ((0.5 * Dstvwo) * ht1 + (Dst - 0.5 * Dstvwo) * ht2) / Dst  # 混合后调节级蒸汽焓值ht
            t = self.a.t_ph(pt * 10, ht)
        else:
            '''部分开阀'''
            pt = Dst * k  # 变工况下调节级后压力
            pst2 = (Dst - 0.75 * Dstvwo) * ((pstb2vwo - ptvwo) / Dstbvwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.a.h_pt(Pst * 10, Ts)  # 主蒸汽焓值
            vst2 = self.a.v_ph(pst2 * 10, hst2)  # 流经部分开阀后的蒸汽的比体积
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht2 = self.a.h_pt(pt * 10, tt2)  # 通过部分开启调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值ht
            # tt = IAPWS_IF97('T_ph', pt, htb - 273.15) #调节级蒸汽温度tt
            '''全开阀'''
            hst1 = self.a.h_pt(Pst * 10, Ts)  # 主蒸汽焓值
            pst1 = 0.95 * Pst  # 蒸汽流经全开调节阀门后压力pst1
            vst1 = self.a.v_ph(pst1 * 10, hst1)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数
            vt1 = (pst1 / pt) ** (1 / nm) * vst1  # 通过全开调节阀蒸汽流经喷嘴动叶后比体积
            tt1 = (1000 * pt * vt1 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过全开调节阀蒸汽流经喷嘴动叶后温度
            ht1 = self.a.h_pt(pt * 10, tt1)  # 通过全开调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值
            '''混合'''
            Dt = Dst  # 调节级后蒸汽流量
            ht = ((0.75 * Dstvwo) * ht1 + (Dst - 0.75 * Dstvwo) * ht2) / Dst  # 混合后调节级蒸汽焓值ht
            t = self.a.t_ph(pt * 10, ht)
        ne = Dst*(hst-ht)
        return Dt, pt, ht, t, ne

    "压力级第一级"
    def pressure_level1(self, Dt, pt, ht,k):
        # p1d = 6.80877688  # 设计工况下第一级抽汽压力
        p1d = 6.90877688 #冬季
        xl = 0.9  # 级间效率
        Dj1d = 583.3333333  # 设计工况下第一级蒸汽流量
        p1 = Dt * k  # 第一级抽汽压力，即第一级出口蒸汽压力
        s1 = self.a.s_ph(pt * 10, ht)  # 级组一入口蒸汽熵
        h11 = self.a.h_ps(p1 * 10, s1)  # 第一级抽汽理想焓值
        h1 = ht - xl * (ht - h11)  # 第一级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps1 = (p1 + Pa) * (1 - EGJCQ / 100) - Pa
        if ps1 > p1:
            D1 = 0
        else:
            D1 = math.sqrt(p1-ps1)*(-0.4972*p1**3+7.227*p1**2-19.44*p1+48.69)  # 第一级抽汽流量
        ne = Dt*(ht-h1)
        return D1, p1, h1, ne

    "压力级第二级"
    def pressure_level2(self, D1, p1, h1, Dt,k):
        # p2d = 4.04  # 设计工况下第二级抽汽压力
        p2d = 4.065 #冬季
        xl = 0.9  # 级间效率
        Dj2d = 550.548348  # 设计工况下第二级蒸汽流量
        Dj2 = Dt - D1  # 第二级蒸汽流量
        p2 = Dj2 * k  # 第二级抽汽压力，即第二级出口蒸汽压力
        s2 = self.a.s_ph(p1 * 10, h1)  # 级组二入口蒸汽熵
        h22 = self.a.h_ps(p2 * 10, s2)  # 第二级抽汽理想焓值
        h2 = h1 - xl * (h1 - h22)  # 第二级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps2 = (p2 + Pa) * (1 - EGJCQ / 100) - Pa
        if ps2 > p2:
            D2 = 0
        else:
            D2 = math.sqrt(p2-ps2)*(-2.75*p2**3+23.6*p2**2-38.75*p2+62.27)
        ne = Dj2*(h1-h2)
        return D2, p2, h2, Dj2, ne

    "压力级第三级"
    def pressure_level3(self, pre, tre, Dre,k):
        # p3d = 2.09  # 设计工况下第三级抽汽压力
        p3d = 2.090 #冬季
        xl = 0.9  # 级间效率
        Dj3d = 520.07  # 设计工况下第三级蒸汽流量
        p3 = Dre * k  # 第三级抽汽压力，即第三级出口蒸汽压力
        hre = self.a.h_pt(pre * 10, tre)  # 级组三入口蒸汽焓值
        s3 = self.a.s_ph(pre * 10, hre)  # 级组三入口蒸汽熵值
        h33 = self.a.h_ps(p3 * 10, s3)  # 第三级抽汽理想焓值
        h3 = hre - xl * (hre - h33)  # 第三级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps3 = (p3 + Pa) * (1 - EGJCQ / 100) - Pa
        if ps3 > p3:
            D3 = 0
        else:
            D3 = math.sqrt(p3-ps3)*(-6.841*p3**3+31.39*p3**2-17.2*p3+22.72)  # 第三级抽汽流量
        Dj3 = Dre - D3  # 第三级蒸汽流量
        ne = Dre*(hre-h3)
        return D3, p3, h3, Dj3, ne, hre

    "压力级第四级"
    def pressure_level4(self, h3, p3, Dj3, k):
        # p4d = 0.99  # 设计工况下第四级抽汽压力
        p4d = 0.916990 #冬季
        xl = 0.9  # 级间效率
        Dj4d = 500.07  # 设计工况下第四级蒸汽流量
        p4 = Dj3 * k  # 第四级抽汽压力，即第四级出口蒸汽压力
        s4 = self.a.s_ph(p3 * 10, h3)  # 级组四入口蒸汽熵
        h44 = self.a.h_ps(p4 * 10, s4)  # 第四级抽汽理想焓值
        h4 = h3 - xl * (h3 - h44)  # 第四级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps4 = (p4 + Pa) * (1 - EGJCQ / 100) - Pa
        if ps4 > p4:
            D4 = 0
        else:
            D4 = math.sqrt(p4-ps4)*(-64.69 *p4**3+144.7*p4**2-47.91*p4+35.01)  # 第四级抽汽流量
        ne = Dj3*(h3-h4)
        Dj4 = Dj3 - D4  # 第四级蒸汽流量

        return D4, p4, h4, Dj4, ne

    "压力级第五级"
    def pressure_level5(self, h4, p4, Dj4, k):
        # p5d = 0.531  # 设计工况下第五级抽汽压力
        p5d = 0.53591 #冬季
        xl = 0.9  # 级间效率
        Dj5d = 485.54  # 设计工况下第五级蒸汽流量
        p5 = Dj4 * k  # 第五级抽汽压力，即第五级出口蒸汽压力
        s5 = self.a.s_ph(p4 * 10, h4)  # 级组五入口蒸汽熵
        h55 = self.a.h_ps(p5 * 10, s5)  # 第五级抽汽理想焓值
        h5 = h4 - xl * (h4 - h55)  # 第五级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps5 = (p5 + Pa) * (1 - EGJCQ / 100) - Pa
        if ps5 > p5:
            D5 = 0
        else:
            D5 = math.sqrt(p5-ps5)*(-365.1*p5**3+529.1*p5**2-68.05 *p5+56.83)  # 第五级抽汽流量
        t5 = self.a.t_ph(p5 * 10, h5)
        ne = Dj4*(h4-h5)
        Dj5 = Dj4 - D5  # 第五级蒸汽流量

        return D5, p5, h5, Dj5, ne

    "压力级第六级"
    def pressure_level6(self, h5, p5, Dj5, k):
        # p6d = 0.2152712  # 设计工况下第六级抽汽压力
        p6d = 0.21902712 #冬季
        xl = 0.9  # 级间效率
        Dj6d = 450.07  # 设计工况下第六级蒸汽流量
        p6 = Dj5 * k  # 第六级抽汽压力，即第六级出口蒸汽压力
        s6 = self.a.s_ph(p5 * 10, h5)  # 级组六入口蒸汽熵
        h66 = self.a.h_ps(p6 * 10, s6)  # 第六级抽汽理想焓值
        h6 = h5 - xl * (h5 - h66)+110  # 第六级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps6 = (p6 + Pa) * (1 - EGJCQ / 100) - Pa
        if ps6 > p6:
            D6 = 0
        else:
            D6 = math.sqrt(p6-ps6)*(9233*p6**3-4349*p6**2+1313*p6+2.72)  # 第六级抽汽流量
        t6 = self.a.t_ph(p6 * 10, h6)
        ne = Dj5*(h5-h6)
        Dj6 = Dj5 - D6  # 第六级蒸汽流量
        return D6, p6, h6, Dj6, ne

    "压力级第七级"
    def pressure_level7(self, h6, p6, Dj6, k):
        # p7d = 0.0864  # 设计工况下第七级抽汽压力
        p7d = 0.081204 #冬季
        xl = 0.9  # 级间效率
        Dj7d = 420.95  # 设计工况下第七级蒸汽流量
        p7 = Dj6 * k  # 第七级抽汽压力，即第七级出口蒸汽压力
        s7 = self.a.s_ph(p6 * 10, h6)  # 级组七入口蒸汽熵
        h77 = self.a.h_ps(p7 * 10, s7)  # 第七级抽汽理想焓值
        h7 = h6 - xl * (h6 - h77)  # 第七级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps7 = (p7 + Pa) * (1 - EGJCQ / 100) - Pa
        if ps7 > p7:
            D7 = 0
        else:
            D7 = math.sqrt(p7-ps7)*(5608*p7**3-604*p7**2+2841*p7-54.62)  # 第七级抽汽流量
        t7 = self.a.t_ph(p7 * 10, h7)
        ne = Dj6*(h6-h7)
        Dj7 = Dj6 - D7  # 第七级蒸汽流量
        return D7, p7, h7, Dj7, ne

    "末级"
    def final_level(self, h7, p7, Dj7, Pc):
        xl = 0.9
        Dc = Dj7  # 末级蒸汽流量
        sm = self.a.s_ph(p7 * 10, h7)  # 末级入口蒸汽熵
        hm = self.a.h_ps(Pc * 10, sm)  # 末级抽汽理想焓值
        hc = h7 - xl * (h7 - hm)  # 末级排汽焓值
        ne = Dc*(h7-hc)
        T = self.a.t_ph(Pc*10, hc)
        Tl_p = self.a.tsat_p(Pc*10)
        return hc, Dc, ne, T, Tl_p

    def simulation(self, Ts, Pst, Dst, P2_real,T2_real, Pc, k):
        Dt, Pt, ht, t, ne0= self.regulating_level(Ts, Pst, Dst,k[0])
        D1j = Dt
        D1, p1, h1, ne1 = self.pressure_level1(Dt, Pt, ht, k[1])

        D2, p2, h2, Dj2, ne2 = self.pressure_level2(D1, p1, h1,Dt, k[2])
        D3, p3, h3, Dj3, ne3, hre = self.pressure_level3(P2_real, T2_real, Dj2, k[3])
        D4, p4, h4, Dj4, ne4 = self.pressure_level4(h3, p3, Dj3, k[4])
        D5, p5, h5, Dj5, ne5 = self.pressure_level5(h4, p4, Dj4, k[5])
        D6, p6, h6, Dj6, ne6 = self.pressure_level6(h5, p5, Dj5, k[6])
        D7, p7, h7, Dj7, ne7 = self.pressure_level7(h6, p6, Dj6, k[7])
        hc, Dc, ne8, T, Tl_p= self.final_level(h7, p7, Dj7, Pc)
        Ne = (ne0+ne1+ne2+ne3+ne4+ne5+ne6+ne7+ne8)*self.ita*0.001
        # return Dt, pt, ht, D1, p1, D1j, h1, D2, p2, h2, D3, p3, h3, Dj3, D4, p4, h4, Dj4,\
        #         D5, p5, h5, Dj5, D6, p6, h6, Dj6, D7, p7, h7, Dj7, hc, Dc, Ne, Tc,Tl_p, delta_h

        return Dt, D1j, Pt, Dc, Ne

"凝汽器"
class Condenser:
    def __init__(self, current_kec):
        self.a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.integrator = Integrator
        self.current_kec = current_kec
        self.Vf_prev = 2.2  # 上一时刻的 Vf，用于风速滤波
        self.Tf = 500.0     # 风速惯性时间常数 [s]
        self.dt = 0.1       # 步长

    def Condenser(self, Pc, n, n_num, Dc, Ta):
        k = 0.04  # 换热系数，W / (m2·℃)
        Vn = 2.5
        Ac = 1440216
        rhoa = 1.043  # 空气密度，kg / m3
        Af = 10731
        Cpa = 1.005  # 空气比热容，kJ / (kg·℃)
        o = 0.9      # 风机集群因子
        Qvd = 550    # 风机体积流量 m³/s
        nd = 110     # 设计转速
        Rs = 461.5   # 气体常数
        Vc = 7700 * 90  # 凝汽器体积，m³

        p = Pc * 10  # MPa → bar
        a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        Hc1 = a.hV_p(p)     # 饱和汽焓
        Tc = a.tsat_p(p)    # 饱和温度
        Hlc = a.hL_p(p)     # 饱和水焓

        # 1. Vf 滤波：一阶惯性
        # Vf_target = Qvd * n_num * o / Af / nd * n
        # self.Vf_prev += (self.dt / self.Tf) * (Vf_target - self.Vf_prev)
        # Vf = self.Vf_prev

        Vf = Qvd * n_num * o / Af / nd * n
        kc = k * (np.power(complex(Vf / Vn), 0.6)).real
        NTU = kc * Ac / rhoa / Vf / Af / Cpa
        Qc = rhoa * Vf * Af * Cpa * (1 - math.exp(-NTU)) * (Tc - Ta)
        Dec = Qc / (Hc1 - Hlc) * self.current_kec
        deta_Pc = Rs * (Tc + 273.15) * (Dc - Dec) / Vc

        return deta_Pc, Dec, Tc, Hc1

    def simulation(self, Pc, n, n_num, Dc, Ta, h):
        deta_Pc, Dec, Tc, Hc1 = self.Condenser(Pc, n, n_num, Dc, Ta)
        Pc = self.integrator.euler(deta_Pc * 0.00001, Pc, h)
        return Pc, Dec, Tc, Hc1

def Net_power(Ne, n_opt, n_num):
    power_consumption_fan = (n_opt / 110) ** 3 * 115 * 0.001 * n_num
    current_net_power = Ne - power_consumption_fan
    return power_consumption_fan, current_net_power

def sim_run(data, k, h, sample_time, current_ita, current_kec, state_result):
    turbine = steam_turbine(current_ita)
    condenser = Condenser(current_kec)

    Pst_DT = data [3]
    Dst_DT = data [4]
    Tst_DT = data[5]
    Pre_DT = data [7]
    Tre_DT = data [8]
    Ta = data [11]
    n = data [12]
    n_num = data [13]

    Pc = state_result[-1] if state_result is not None and len(state_result) > 0 else 0.007

    Ds, D1j, Pt, Dc, Ne = turbine.simulation(Tst_DT, Pst_DT, Dst_DT, Pre_DT, Tre_DT, Pc, k)
    for j in range(int(sample_time / h)):
        Pc, Dec, Tc, Hc1 =  condenser.simulation(Pc, n, n_num, Dc, Ta, h)
    state_result = np.array([Ds, D1j, Pt, Dc, Ne, Pc])
    return state_result

def scope(x, y, title, y_true=None):
    # 设置字体配置
    size = 16
    config = {
        "font.family": "serif",
        "font.size": size,
        "axes.labelsize": size,
        "axes.titlesize": size,
        "figure.titlesize": size,
        "xtick.labelsize": size,
        "ytick.labelsize": size,
        "legend.fontsize": size,
        "legend.title_fontsize": size,
        "mathtext.fontset": "stix",
        "font.serif": ["Times New Roman"],
        "axes.unicode_minus": False,
    }
    rcParams.update(config)

    # 创建图表
    plt.figure(figsize=(10, 4))  # 设置图表大小

    # 绘制曲线
    plt.plot(x, y, color="g", linestyle="-.", label="y")
    if y_true is not None:
        plt.plot(x, y_true, color="b", linestyle="-.", label="y_true")

    # 设置图表属性
    plt.xlabel("time (minute)")
    plt.ylabel("Y")
    plt.title(title)
    plt.legend(loc="best")  # 使用更具描述性的参数值

    # 自动调整布局
    plt.tight_layout()

    # 保存图表并显示
    plt.savefig(f"{title}.png", dpi=300, bbox_inches="tight")  # 提高保存图片质量
    plt.show(block=True)

    return

def De_cal(p1_real,p2_real,p3_real,p4_real,p5_real,p6_real,p7_real):
    Pa = 0.101
    EGJCQ = 3
    ps1_real = (p1_real + Pa) * (1 - EGJCQ / 100) - Pa
    D1_real = math.sqrt(p1_real-ps1_real)*(-0.4972*p1_real**3+7.227*p1_real**2-19.44*p1_real+48.69)

    ps2_real = (p2_real + Pa) * (1 - EGJCQ / 100) - Pa
    D2_real = math.sqrt(p2_real-ps2_real)*(-0.4972*p2_real**3+7.227*p2_real**2-19.44*p2_real+48.69)

    ps3_real = (p3_real + Pa) * (1 - EGJCQ / 100) - Pa
    D3_real = math.sqrt(p3_real-ps3_real)*(-0.4972*p3_real**3+7.227*p3_real**2-19.44*p3_real+48.69)

    ps4_real = (p4_real + Pa) * (1 - EGJCQ / 100) - Pa
    D4_real = math.sqrt(p4_real-ps4_real)*(-0.4972*p4_real**3+7.227*p4_real**2-19.44*p4_real+48.69)

    ps5_real = (p5_real + Pa) * (1 - EGJCQ / 100) - Pa
    D5_real = math.sqrt(p5_real-ps5_real)*(-0.4972*p5_real**3+7.227*p5_real**2-19.44*p5_real+48.69)

    ps6_real = (p6_real + Pa) * (1 - EGJCQ / 100) - Pa
    D6_real = math.sqrt(p6_real-ps6_real)*(-0.4972*p6_real**3+7.227*p6_real**2-19.44*p6_real+48.69)

    ps7_real = (p7_real + Pa) * (1 - EGJCQ / 100) - Pa
    D7_real = math.sqrt(p7_real-ps7_real)*(-0.4972*p7_real**3+7.227*p7_real**2-19.44*p7_real+48.69)
    return D1_real,D2_real,D3_real,D4_real,D5_real,D6_real,D7_real

def optimize_for_speed(data1, m1, h, sample_time, tolerance, current_ita, current_kec, K, k_coal, k_Dt, state_result=None):
    """
    优化风机转速，根据仿真负荷与实际负荷的差值动态调整主汽流量

    参数:
    max_iterations: 负荷-主汽流量迭代的最大次数
    """
    # k_coal = 0.9574  # 负荷与主汽流量关系参数，需从外部输入
    # b = -20.349  # 负荷与主汽流量关系参数，需从外部输入
    # k_Dt = 0.6039  # 主汽流量与给煤量关系参数，需从外部输入
    # d = 29.2676  # 主汽流量与给煤量关系参数，需从外部输入

    # 输入煤价和电价
    # load_rate = data1[1]/600  # 实际负荷
    # coal_price = 500 if load_rate <= 0.5 else 550 if load_rate <= 0.75 else 625
    coal_price = 550.0/1000  # 煤价 (元/吨)
    electricity_price = 0.35291  # 电价 (元/kWh)
    max_iterations = 50

    # 保存原始数据
    actual_load = data1[1]  # 实际负荷
    actual_back_pressure = data1[2]   # kPa，实际背压
    actual_D = data1[4] # 原始主汽流量
    actual_fan_speed = data1[12]  # 实际风机转速

    data1[12] = m1  # 设置风机转速
    n_num = data1[13]  # 风机运行台数

    # 主汽流量-负荷迭代调整
    for iteration in range(max_iterations):
        # 运行仿真计算当前背压和负荷
        state_result = sim_run(data1, K, h, sample_time, current_ita, current_kec, state_result=state_result)
        current_Pc = state_result[-1]*1000  # 当前背压 (MPa-→kPa)
        current_Ne = state_result[4]  # 当前负荷 (MW)

        delta_Ne = data1[1] - current_Ne  # 仿真与实际负荷差
        delta_Ne = np.clip(delta_Ne, -20, 20)
        if delta_Ne > 20:
            delta_Ne = 20
        if delta_Ne < -20:
            delta_Ne = -20

        dD = (k_Dt * delta_Ne)

        # 根据负荷差值调整主汽流量
        D_adj = data1[4] + dD

        # 更新数据中的主汽流量
        data1[4] = D_adj

        # 判断是否收敛
        if abs(delta_Ne / actual_load) < tolerance or iteration == max_iterations - 1:
            # print(f"风机转速 {m1}: 负荷={current_Ne:.2f} MW, 背压={current_Pc:.2f} kPa, "
            #       f"主汽流量={D_adj:.2f} t/h, 负荷差值={load_diff:.2f} MW")# 打印迭代信息
            break

    # 最终背压和负荷
    final_back_pressure = current_Pc
    final_load = current_Ne

    # 基于主汽流量变化计算给煤量减少量，给煤量与主汽流量成线性关系: Coal = c * D + d
    coal_reduction = k_coal * (actual_D - D_adj)  # 给煤量减少量

    # 计算给煤节省的成本 (元/分钟 或 元)
    coal_cost_saving = coal_reduction * 60 * coal_price

    # 计算净功率
    power_consumption_fan, net_power = Net_power(float(final_load), m1, n_num)

    # 计算风机增加的功耗成本 (元/分钟)
    power_cost_increase = (power_consumption_fan - data1[15]) *1000 / 60* electricity_price   # 转换为 kW

    # 经济效益计算
    objective_value = coal_cost_saving - power_cost_increase

    return [
        round(data1[11], 2),  # 环境温度
        round(net_power, 2),  # 最大净功率
        round(data1[16], 2),  # 实际净功率
        round(final_load, 2),  # 优化负荷
        round(actual_load, 2),  # 实际负荷
        round(final_back_pressure, 2),  # 最优背压
        round(actual_back_pressure, 2),  # 实际背压
        round(m1, 2),  # 优化风机转速
        round(actual_fan_speed, 2),  # 实际风机转速
        round(data1[13], 2),  # 风机运行台数
        round(power_consumption_fan, 2),  # 优化风机功耗
        round(data1[15], 2),  # 实际风机功耗
        round(D_adj, 2),  # 调整后的主汽流量
        round(actual_D, 2),  # 原始主汽流量
        round(coal_reduction, 2),  # 给煤量减少量
        round(coal_cost_saving, 2),  # 给煤节省成本 (元/分钟)
        round(power_cost_increase, 2),  # 风机增加成本 (元/分钟)
        round(objective_value, 2)  # 目标函数值 (元/分钟)
    ], state_result

def backpressure_optimize(data, output_filepath, sample_time=60, h=0.1,
                      tolerance=0.025, speed_change_limit=16.5, speed_step=1):
    """
    单线程背压优化函数，保存单个样本在所有风机转速下的运行结果，并取出最优结果

    参数:
    speed_change_limit: 相邻样本风机转速最大变化量
    speed_step: 遍历风机转速时的步长
    """
    result_list = []  # 存储所有样本的最优结果
    all_samples_results = []  # 存储所有样本在所有风机转速下的完整结果

    prev_optimal_speed = None  # 上一时刻的最优风机转速

    # 初始化模型参数
    current_ita = 0.94  # 汽轮机效率
    current_kec = 1.5  # 流量系数
    ita_ = [current_ita]
    kec_ = [current_kec]

    # 初始化状态结果
    state_result = [data[0,4], data[0,4], data[0,18], data[0,10], data[0,1], data[0,2]/1000]

    # 遍历所有样本
    for i in tqdm(range(len(data))):
        if i == 0:
            D = np.asarray([[100, 100, 100, 100, 100, 100, 100, 100]])
            P = np.asarray([[100, 100, 100, 100, 100, 100, 100, 100]])
            K = [12.3 / 513.6, 6.9 / 583, 4 / 550, 2 / 520, 0.9 / 500, 0.53 / 485, 0.21 / 450, 0.08 / 420]

        if np.any(data[i,4:12]< 0.001) :
            continue

        actual_load = data[i, 1]  # 当前样本的实际负荷
        actual_back_pressure = data[i, 2]  # 当前样本的实际背压 (kPa)

        # ---- 1. 参数更新 ----
        # 设置收敛条件
        back_pressure_tolerance = 0.01  # 背压误差容忍度 (1%)
        max_iterations = 50  # 最大迭代次数
        param_converged = False
        iteration = 0

        while not param_converged and iteration < max_iterations:
            state_result = sim_run(data[i], K, h, sample_time, current_ita, current_kec, state_result=state_result)

            result_Pc = state_result[-1] * 1000  # 背压 (kPa)
            result_Ne = state_result[4]  # 负荷 (MW)

            # 计算模型误差
            error_Ne = result_Ne - actual_load
            error_Pc = result_Pc - actual_back_pressure

            # 自适应更新模型参数
            gamma1 = 0.05  # 背压参数更新步长
            current_kec += gamma1 * error_Pc
            current_kec = np.clip(current_kec, 0.5, 5)
            kec_.append(current_kec)

            gamma2 = -0.0005  # 负荷参数更新步长
            current_ita += gamma2 * error_Ne
            current_ita = np.clip(current_ita, 0.75, 1.0)
            ita_.append(current_ita)

            # 检查收敛条件
            if abs(error_Pc)/actual_back_pressure < back_pressure_tolerance:
                param_converged = True

            iteration += 1

        Ds, D1j, Pt, Dc, Ne, Pc = (state_result[item] for item in range(len(state_result)))
        D1_real, D2_real, D3_real, D4_real, D5_real, D6_real, D7_real = (De_cal(data[i, 19], data[i, 20],
                                                                        data[i, 21], data[i, 22], data[i, 23], data[i, 24],data[i, 25])) #抽汽压力
        P_save = np.asarray([data[i, 18], data[i, 19], data[i, 20], data[i, 21], data[i, 22], data[i, 23], data[i, 24], data[i, 25]])
        Dj2_real = D1j - D1_real
        Dj3_real = Dj2_real - D2_real
        Dj4_real = Dj3_real - D3_real
        Dj5_real = Dj4_real - D4_real
        Dj6_real = Dj5_real - D5_real
        Dj7_real = Dj6_real - D6_real

        buf = 10

        D_save = np.asarray([[Ds, D1j, Dj2_real, Dj3_real, Dj4_real, Dj5_real, Dj6_real, Dj7_real]])
        if D.shape [0] <= buf:
            D = np.vstack((D, D_save))
            P = np.vstack((P, P_save))
        if D.shape [0] > buf:
            D = D [1:, :]
            P = P [1:, :]
            k0 = turbine_coff(D [:, 0].reshape(1, -1), P [:, 0].reshape(1, -1))
            k1 = turbine_coff(D [:, 1].reshape(1, -1), P [:, 1].reshape(1, -1))
            k2 = turbine_coff(D [:, 2].reshape(1, -1), P [:, 2].reshape(1, -1))
            k3 = turbine_coff(D [:, 3].reshape(1, -1), P [:, 3].reshape(1, -1))
            k4 = turbine_coff(D [:, 4].reshape(1, -1), P [:, 4].reshape(1, -1))
            k5 = turbine_coff(D [:, 5].reshape(1, -1), P [:, 5].reshape(1, -1))
            k6 = turbine_coff(D [:, 6].reshape(1, -1), P [:, 6].reshape(1, -1))
            k7 = turbine_coff(D [:, 7].reshape(1, -1), P [:, 7].reshape(1, -1))
            K = np.asarray([k0, k1, k2, k3, k4, k5, k6, k7]).reshape(-1)

        if i > 360:
            k_coal = (data[i-360:i, -1].reshape(1, -1) @ np.linalg.pinv(
                np.concatenate((data[i-360:i, 4].reshape(1, -1), np.ones((1, 360))))))  # .reshape(-1)
            k_Dt = (data[i-360:i, 4].reshape(1, -1) @ np.linalg.pinv(
                np.concatenate((data[i-360:i, 1].reshape(1, -1), np.ones((1, 360))))))  # .reshape(-1)
            k_coal = k_coal[0, 0]
            k_Dt = k_Dt[0, 0]
        else:
            k_coal = 0.1929
            k_Dt = 0.8624

        # ---- 根据负荷率计算背压上下限 ----
        full_load = 600
        load_rate = actual_load / full_load  # 当前数据点的负荷率

        max_back_pressure = 41.8182 * load_rate + 16.6364 if 0.2 < load_rate <= 0.75 else 25 if load_rate <= 0.2 else 48
        min_back_pressure = 8.75 * load_rate - 0.75 if load_rate > 0.1 else 0

        # ---- 2. 背压优化 ----
        start_t = time.time()

        # min_speed, max_speed = 44, 121

        # 确定风机转速搜索范围
        if prev_optimal_speed is None:
            # 第一个样本使用完整范围
            min_speed, max_speed = 44, 120.5
            # min_speed, max_speed = 70, 100
        else:
            # 后续样本考虑相邻时刻变化率约束
            min_speed = max(44, prev_optimal_speed - speed_change_limit)
            max_speed = min(120.5, prev_optimal_speed + speed_change_limit)

        print(f"样本 {i + 1}/{len(data)}: 风机转速搜索范围 [{min_speed:.1f}, {max_speed:.1f}] rpm")

        # 生成要遍历的风机转速列表
        speed_range = np.arange(min_speed, max_speed + speed_step, speed_step)
        print(f"需要评估的风机转速点数量: {len(speed_range)}")

        # 存储当前样本在所有风机转速下的结果
        current_sample_results = []
        best_result = None
        best_objective = -float('inf')
        best_speed = None

        # 遍历每个风机转速，评估经济效益
        prev_valid_bp = float('inf')  # 初始化上一个有效背压为无穷大
        for m1 in speed_range:
            data_copy = data[i, :].copy()

            if m1 == min_speed:
                result, state_result = optimize_for_speed(data_copy, m1, h, sample_time, tolerance, current_ita,
                                                           current_kec, K, k_coal, k_Dt, state_result=state_result)
            else:
                result, state_result = optimize_for_speed(data_copy, m1, h, sample_time, tolerance, current_ita,
                                                           current_kec, K, k_coal, k_Dt, state_result=state_result)
            if m1 > min_speed:
                bp = result[5]  # 最优背压 (kPa)
                objective = result[-1]  # 目标函数值 (元/分钟)

                # 添加风机转速和样本索引到结果
                extended_result = [i + 1] + result  # 添加样本编号和风机转速
                current_sample_results.append(extended_result)

                # 检查背压是否在合理范围内
                if min_back_pressure <= bp <= max_back_pressure:
                    # 检查背压是否不再下降
                    # if bp >= prev_valid_bp:
                    #     print(f"风机转速为{m1:.2f} 时，背压{bp:.2f} kPa不再下降（上一个有效背压: {prev_valid_bp:.2f} kPa）")
                    #     break  # 背压不再下降时跳出循环

                    # prev_valid_bp = bp  # 更新上一个有效背压

                    if objective > best_objective:
                        best_objective = objective
                        best_result = extended_result
                        best_speed = m1
                else:
                    # print(f"风机转速为{m1:.2f} 时，背压{bp:.2f} kPa超出合理范围")
                    continue  # 背压超出范围时跳出循环

            # 打印进度
            if (m1 - min_speed) % (speed_step * 10) == 0:
                print(f"  已评估 {int((m1 - min_speed) / speed_step + 1)}/{len(speed_range)} 个转速点")

        # 将当前样本的所有结果添加到总结果列表
        all_samples_results.extend(current_sample_results)

        end_t = time.time()
        print(f'样本 {i + 1}/{len(data)} 优化完成，耗时{(end_t - start_t):.2f}秒，'
              f'最优风机转速: {best_speed} rpm，经济效益: {best_objective:.2f} 元/分钟')

        # 记录最优结果
        if best_result is not None:
            print(f"已找到最优解（数据点 {i + 1}/{len(data)}）")

            result_list.append(best_result)
            prev_optimal_speed = best_speed  # 保存最优风机转速

        if i%360==0 or i==len(data)-1:

            # 保存所有样本在所有风机转速下的完整结果
            all_results_df = pd.DataFrame(all_samples_results, columns=[
                '样本编号',
                '环境温度',
                '最大净功率',
                '实际净功率',
                '优化负荷',
                '实际负荷',
                '最优背压',
                '实际背压',
                '优化风机转速',
                '实际风机转速',
                '风机运行台数',
                '优化风机功耗',
                '实际风机功耗',
                '调整后主汽流量',
                '原始主汽流量',
                '给煤量减少量',
                '给煤节省成本 (元/分钟)',
                '风机增加成本 (元/分钟)',
                '经济效益 (元/分钟)',
            ])

            all_results_file = output_filepath.replace('.xlsx', '_所有结果.xlsx')
            try:
                all_results_df.to_excel(all_results_file, index=False)
                print(f"所有样本在所有风机转速下的完整结果已保存至：{all_results_file}")
            except Exception as e:
                print(f"保存完整结果时出错: {e}")

            # 绘制自适应系数变化图
            plt.figure(figsize=(12, 7))

            plt.subplot(2, 1, 1)
            plt.plot(ita_, 'r-')
            plt.ylabel('汽轮机效率')
            plt.title('汽轮机效率')

            plt.subplot(2, 1, 2)
            plt.plot(kec_, 'g-')
            plt.ylabel('流量系数')
            plt.title('流量系数变化趋势')

            plt.tight_layout()
            plt.savefig('图片/自适应系数.png', dpi=300)
            # plt.show()

            # 将最优结果保存到一个DataFrame中
            result_df1 = pd.DataFrame(result_list, columns=[
                '样本编号',
                '环境温度',
                '最大净功率',
                '实际净功率',
                '优化负荷',
                '实际负荷',
                '最优背压',
                '实际背压',
                '优化风机转速',
                '实际风机转速',
                '风机运行台数',
                '优化风机功耗',
                '实际风机功耗',
                '调整后主汽流量',
                '原始主汽流量',
                '给煤量减少量',
                '给煤节省成本 (元/分钟)',
                '风机增加成本 (元/分钟)',
                '经济效益 (元/分钟)',
            ])

            print("优化结果维度为：", len(result_list))

            try:
                result_df1.to_excel(output_filepath, index=False)
                print(f"最优结果已保存至：{output_filepath}")
            except Exception as e:
                raise IOError(f"保存 Excel 文件时出错: {e}")

    return result_list

if __name__ == '__main__':
    ini_num = 1440*7+12156
    # sim_num = data.shape[0]
    sim_num = 2

    # df = pd.read_excel(r'H:\2023数据\8.xlsx')
    # data = df.to_numpy()
    # file_path = r'全部数据\1分钟\10月_处理后.npy'
    # data = np.load(file_path)
    data = np.load("D:/桌面文件/桌面文件/太原理工大学/毕业设计/背压优化及控制/数据/8月/npy文件/8月数据.npy",
                   allow_pickle=True)
    output_filepath = "D:/桌面文件/桌面文件/太原理工大学/毕业设计/背压优化及控制/数据/8月/excel数据/冷端背压优化终版结果8.8-5.31.xlsx"
    data_sim = data[ini_num:ini_num+sim_num, :]
    # data_sim = data_sim.astype(float)
    result_list = backpressure_optimize(data_sim, output_filepath, sample_time=60, h=0.1, tolerance=0.015)
